args <- commandArgs(TRUE)
setwd(args[1])
library(survival)
#reading the input data(here genotype)
geno<-read.table("temp_input",sep=" ",head=T)
#geno<-read.table("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/input_withcorrecthead.tped",sep=" ",head=T)
rownames(geno)<-geno[,2]

#reorder samples according Martingae method
ordersamples<-read.table("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/ordered_sampes_Martingale.txt")
ordersamples<-as.vector(ordersamples[,1])
ordersamples<-gsub("-",".",ordersamples)
ordersamples<-c("chr","rsid","dist","pos",ordersamples)
geno<-geno[,ordersamples]

#reading the clinical data
clin<-read.delim("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/clinical" ,sep="\t",quote="\"",dec=".",fill=TRUE,comment.char="",
header=TRUE,
col.names=c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE"),as.is=c(TRUE,TRUE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE),colClasses=c("character","numeric","factor","numeric","numeric","numeric","integer","integer","integer","factor"))
clin$SUBJECT<-gsub("-",".",clin$SUBJECT)
samplesize<-20
num<-as.integer((length(geno[1,])-4)/samplesize)
num<-20
AUCrank<-matrix(nrow=nrow(geno), ncol=num)
rownames(AUCrank)<-geno[,2]
#colgeno<-colnames(geno)
#colnames(hazardrank)<-colgeno[5:length(colgeno)]
#predgmes(dt) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","EXPR")ene<-c()
i<-1

#looping through each snp
while(i<nrow(geno)+1)
#while(i<3)
{
	geno1<-geno[i,5:ncol(geno)]
#generating 20 different order vectors based on the input vector(current snp in the loop)
#as d change sub iterates through all the samples
	d<-1
	while(d<num+1)
	{			
		s<-d;
		sub<-c()
		main<-c()
		if(d > 1)
		{
			main<-seq(1,d-1,1)
		}
		nu<-1
		while(s<length(geno1)+1)
		{
			sub[nu]<-s
			sequ<-c()		
			if(s+1 <= length(geno1))
			{	
				sequ<-seq(s+1,s+num-1,1)
				if(s+num-1 > length(geno1))
				{
					sequ<-seq(s+1,length(geno1),1)
				}
			}
			main<-c(main,sequ)
			nu<-nu+1	
			s<-s+num
		}
#main,sub contains (main 19 parts) (sub(1 part)

		#print(sub)
#creating survival model using expr and survival time of the main vector(19 parts out of 20) samples
		genomain<-geno1[main]
		genosub<-geno1[sub]

		dt <- merge(clin,t(genomain),by.y="row.names",by.x="SUBJECT")
		colnames(dt) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","EXPR")
		fit1 <- coxph(Surv(OS,ALIVEOS) ~ EXPR+AGE,data=dt)
		#pred<-geno1
		#pred <- predict(fit1, data.frame(cbind(genotype=geno1)), type="response")

#with the above model predicting the  values for all the samples (all 20 parts)
		dt1 <-merge(clin,t(geno1),by.y="row.names",by.x="SUBJECT")
		colnames(dt1) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","EXPR")
		pred <- predict(fit1, dt1,type="lp")
		#order_rank<-order(pred)
                #sorted_ranks<-seq(1,length(order_rank),1)
                #order_of_index<-order(order_rank)
                #anted_ranks<-sorted_ranks[order_of_index]
#calculating hazard ranks calculating AUC
		wanted_ranks<-rank(pred)
		wr1<-matrix(nrow=1, ncol=length(wanted_ranks))
		colnames(wr1)<-colnames(geno1)
		rownames(wr1)<-rownames(geno1)
		wr1[1,]<-wanted_ranks
		dt <- merge(clin,t(wr1),by.y="row.names",by.x="SUBJECT")
		colnames(dt) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","RANK")
        	#surv<-survConcordance(Surv(OS,ALIVEOS)~RANK,data=dt)
		#AUC<-surv$concordance
		#Hugues code start
		dt1<-dt[rev(sort.list(dt[,11])), ]
		i1<-1
		hazard_rank<-dt1[,11]
		N<-nrow(dt1)
		survival_time<-dt1[,4]
		k<-0
		SensArr<-c()
		oneMSpecArr<-c()
		nsens<-1
		SensArr[nsens]<-0
		oneMSpecArr[nsens]<-0
		nsens<-nsens+1
		while(i1<nrow(dt1)+1)
		{
			if(k != dt1[i1,11])
			{
				#pi<-i
				t<-dt1[i1,4]
				#coln<-colnames(geno)
				c<-hazard_rank[i1]
				Pr_M_ge_c<-length(which(hazard_rank>=c))/N
				Pr_T_le_t<-length(which(survival_time<=t))/N	
				Pr_T_le_t_given_M_ge_c<-length(which(survival_time<=t & hazard_rank>=c))/N
				Pr_T_gt_t_given_M_lt_c<-length(which(survival_time>t & hazard_rank<c))/N
				Sens<-(Pr_T_le_t_given_M_ge_c*Pr_M_ge_c)/Pr_T_le_t
				if(Pr_T_le_t == 1) {
					Spec<-0
				} else
				 {
				        Spec<-Pr_T_gt_t_given_M_lt_c*(1-Pr_M_ge_c)/(1-Pr_T_le_t)
				}
				oneMSpec<-1-Spec
				SensArr[nsens]<-Sens
				oneMSpecArr[nsens]<-oneMSpec
				nsens<-nsens+1	

				#print(c("Y",y,"dx",dx))
				k<-dt1[i1,11]
			}
			i1<-i1+1
		}
		SensArr[nsens]<-1
		oneMSpecArr[nsens]<-1
		#AUC<-0
		#prev_sens<-0
		#prev_oneMSpec<-0
		#sorted_sens<-order(oneMSpecArr)
		#i1<-1
		#while(i1<=nsens) 
		#{	
		#	Sens=SensArr[sorted_sens[i1]]
		#	oneMSpec=oneMSpecArr[sorted_sens[i1]]
		#	dx<-oneMSpec-prev_oneMSpec
		#	y<-(Sens+prev_sens)/2
		#	AUC<-AUC+y*dx
		#	#print(c("Y",y,"dx",dx))
		#	prev_oneMSpec<-oneMSpec
		#	prev_sens<-Sens
		#	i1<-i1+1
		#}
		#print(c("AUC",AUC))
		trap.rule <- function(x,y) sum(diff(x)*(y[-1]+y[-length(y)]))/2
		sorted_sens<-order(oneMSpecArr)
		AUC<-trap.rule(SensArr[sorted_sens],oneMSpecArr[sorted_sens])

		
		#Hugues code end
		AUCrank[i,d]<-AUC
		#hazardrank[i,sub]<-wanted_ranks[sub]
		d<-d+1
	}
	#print(hazardrank[i,])
	print(i)
	i<-i+1
}

#colnames(hazardrank)<-colnames(geno)
#write.table(hazardrank,"hazardrank.txt",quote=FALSE,sep="\t",col.names=TRUE,row.names=TRUE)
#tnum<-1
#AUCrank<-matrix(nrow=nrow(hazardrank), ncol=1)
#colnames(AUCrank)<-c("AUC")
#rownames(AUCrank)<-rownames(hazardrank)
#while(tnum<nrow(hazardrank)+1)
#{
 #       dt <- merge(clin,hazardrank[tnum,],by.y="row.names",by.x="SUBJECT")
#        colnames(dt) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","RANK")
 #       surv<-survConcordance(Surv(OS,ALIVEOS)~RANK,data=dt)
#        AUC<-surv$concordance
#        #print(c("AUC",AUC))
#        AUCrank[tnum,1]<-AUC
#	print(c(tnum,AUCrank[tnum,1]))
#        tnum =tnum+1
#}
write.table(AUCrank,"AUC.txt",quote=FALSE,sep="\t",col.names=TRUE,row.names=TRUE)

